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Abstract 

We consider the initial-value problem for a system of coupled Boussinesq equations on the 
infinite line for localised or sufficiently rapidly decaying initial data, generating sufficiently 
rapidly decaying right- and left-propagating waves. We study the dynamics of weakly nonlinear 
waves, and using asymptotic multiple-scales expansions and averaging with respect to the fast 
time, we obtain a hierarchy of asymptotically exact coupled and uncoupled Ostrovsky equations 
for unidirectional waves. We then construct a weakly nonlinear solution of the initial-value 
problem in terms of solutions of the derived Ostrovsky equations within the accuracy of the 
governing equations, and show that there are no secular terms. When coupling parameters are 
equal to zero, our results yield a weakly nonlinear solution of the initial-value problem for the 
Boussinesq equation in terms of solutions of the initial- value problems for two Korteweg-de Vries 
equations, integrable by the Inverse Scattering Transform. We also perform relevant numerical 
simulations of the original unapproximated system of Boussinesq equations to illustrate the 
difference in the behaviour of its solutions for different asymptotic regimes. 

Keywords: Coupled Boussinesq equations; Ostrovsky equation; Asymptotic multiple-scales 
expansions; Averaging; Initial-value problem 

1 Introduction 

The Ostrovsky equation 

(r] t + ur\r] x + jrj xxx ) x = Xrj 

is a modification of the Korteweg-de Vries (KdV) equation for the study of oceanic waves, which 
takes into account the effect of background rotation p]. It is well known that the rotation in 
the oceanographic problem (7A > 0) eliminates the solitary wave solutions of the Korteweg- 
de Vries equation through the terminal radiation damping [2 [3] . The numerical simulations in 
Refs. |H IS] have shown that a localised wave packet emerges as a stable dominant solution of the 
Ostrovsky equation. In an independent study [6], it was established that stable envelope solitons 
play a central role in the dynamics of a modified Toda lattice with an additional linear term, 
which can be related to the two-directional generalisation of the Ostrovsky equation derived 



*Accepted for publication in Wave Motion (2011) doi: 10. 1016/j.wavemoti. 2011. 04.003 

t Corresponding author. Tel: +44 (0)1509 228202. Fax: +44 (0)1509 223969. 



in Ref. [7]. The weakly nonlinear description of the emerging wave packet for the Ostrovsky 
equation in terms of a higher-order nonlinear Schrodinger equation has been developed in Ref. 
[5], linking the wavenumber of the carrier wave with the extremum of the group velocity [6]. 

In this paper, we are concerned with the construction of a weakly nonlinear solution of the 
initial-value problem for a system of coupled regularised Boussinesq (cRB) equations [8] 

Utt ~ Uxx = U X U XX + Uttxx -S(U-W), 

w tt - c 2 w xx = aw x w xx + I3w ttxx + l{u - w). (1) 

The regularised version of the Boussinesq equation is preferable from the viewpoint of numerical 
simulations due to suppression of the short wave instability (see Refs. [HI HO]), although this 
version is not integrable by the Inverse Scattering Transform [T2] . From the viewpoint of our 
developed analytical approach, this is not essential, and we could have worked with any version 
of coupled Boussinesq equations. Also, within our approach, generalisation of the derivations to 
the case of three and more equations of this type is straightforward, and we do not discuss it in 
this paper, although the detailed study of various physical effects is interesting. In the context 
of waves in solids, Boussinesq-type equations have been derived, for example, for nonlinear 
waves in solid waveguides [131 CES] an d for waves in microstructured solids [HI [T7] (further 
references can be found in Ref. [TH] and Refs. [HI [23 

System has been recently derived as an accurate asymptotic model for long nonlinear 
longitudinal waves in a layered waveguide with a soft bonding layer using a complex nonlinear 
layered lattice model [8j. The system may be also obtained as a continuum limit for a model of 
two one-dimensional Fermi-Pasta-Ulam (FPU) chains [22] with weak coupling between them. 
System ([I]) is Lagrangian with the Lagrangian density 
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It has three conservation laws A\ + B l x = 0, % = 1, 2, 3, with the densities 
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A =u t + -w t , 
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,3 8 P$ 
A = u t u x + -w t w x + u tx u xx H w tx w xx 

7 7 

(for details see Ref. [8]). 

Differentiating ([!]) with respect to x, and denoting u x = /, w x = g, we obtain 

ftt fxx = 2^(/ )xx ~\~ fttxx ${f d): 

gtt - c 2 g xx = ^a(g 2 ) xx + (3g ttxx + j(f - g) (2) 

(uncoupled equations in this form are sometimes called "regularised long wave equations" and 
"improved bad Boussinesq equations" ) . In what follows we will consider solutions of this system 
instead of the solutions of the original system ([!]). We refer to both systems ([T]) and ^ as 
coupled regularised Boussinesq (cRB) equations, since system ^ is obtained by differentiation 
of system ([l]) . 

We are interested in constructing a weakly nonlinear solution of the initial-value problem 
(IVP) for system @. We use asymptotic multiple-scales expansions of the type used in the 
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study of oblique interaction of solitary waves in Refs. \23\ [24"] (see also references therein and 
in Ref. [25j ) . Recently, we developed a scheme based on this type of asymptotic expansions, 
which allowed us to solve a weakly nonlinear wave scattering problem |26j , formulated in terms 
of a Boussinesq-type equation with piecewise-constant coefficients subject to two continuity 
conditions across the jump and some natural radiation conditions. In this paper, we first use 
the procedure of averaging with respect to the fast time to obtain a hierarchy of asymptotically 
exact coupled and uncoupled Ostrovsky equations for the cases when the characteristic linear 
speeds of the two wave operators in Eqs. ^ are close or essentially different (strong or weak 
interactions in the terminology of Refs. |23 | 124 1 [25]). More precisely, to leading order we derive 
four uncoupled Ostrovsky equations when c — 1 = 0(1), but two coupled systems of Ostrovsky 
equations when c — 1 = 0(e). Then, we show how to construct the weakly nonlinear solution 
of the IVP in terms of solutions of the derived Ostrovsky equations within the accuracy of the 
Eqs. ([2]). We also establish that corrections to the leading-order terms are nonsecular due to a 
special property of solutions of the Ostrovsky equation. In the absence of coupling (5 = 7 = 0), 
these results yield a weakly nonlinear solution of the initial-value problem for the Boussinesq 
equation in terms of solutions of the initial-value problems for two Korteweg-de Vries (KdV) 
equations. Finally, we perform numerical simulations of the original unapproximated system 
([2]) to show the difference in the asymptotic behaviour of its solutions, when initial conditions 
are taken in the form of co-propagating solitary waves of the uncoupled equations. The results 
support our theory. 

The word 'hierarchy' is used here to reflect on the growing complexity of the leading order 
asymptotic models. In particular, generalisation to the case of TV" coupled Boussinesq-type 
equations (describing, for example, long longitudinal waves in TV-layered elastic waveguide, or 
waves in N coupled FPU chains), with N characteristic speeds close to each other, will lead to 
a system of N coupled Ostrovsky equations. 



2 Dispersion curve and solitary waves 

In the symmetric case, when c = a = (3 = 1, system ([2]) admits a reduction g = f, where / 
satisfies the equation 

ftt - f xx = -(f 2 U + ftt xx . (3) 
Eq. ([3]) has particular solutions in the form of solitary waves: 

,x — vt 



f = A sech 2 - 



A 



where A = 3(v 2 — 1), A = —^=. However, these pure or classical solitary wave solutions, 
rapidly decaying to zero in their tail regions, are structurally unstable and are replaced with 
radiating solitary waves [5], i.e. a solitary wave radiating a co- propagating one-sided oscillatory 
tail, using the terminology of Refs. ES] • There have been extensive studies of generalised 

and radiating solitary waves, especially in the context of fluid mechanics (e.g., [30] - [36j), but 
the models were different from Eqs. Q. In particular, long- wave ripples are radiated by solitons 
in Eqs. Q, due to the type of coupling terms in the equations, and the resulting structure of 
the dispersion relation. 

We consider system ([2]), and assume that coefficients are perturbed compared to the sym- 
metric case above, but remain positive. The linear dispersion relation has the form 

[k 2 (l - p 2 ) - k 4 p 2 + 5} [k 2 (c 2 -p 2 )- f3k A p 2 + l] = 1$, 

where k is the wavenumber and p is the phase speed, and was analysed in Ref. [8]. A typical 
linear dispersion curve of Eqs. ([2]) is shown in Fig(T] Significant difference with the linear 
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dispersion curve of the reduction ([3]) consists in the appearance of the second (upper) branch, 
going to infinity as k — > 0, and approaching zero, remaining above the lower branch, as k — > oo. 




k k 

Figure 1: (a) Two branches of the linear dispersion curve of Eqs. gforc = 1.05, = 1, 5 = 7 = 0.01 
and intersection with p = 1.3 (horizontal line) and (b) Linear dispersion curve of the reduction ^ in the 
symmetric case c = l,/3 = l,<5 = 7 = 0.01. 



The classical or pure solitary waves of the single Boussinesq equation ([3]) arise as a bifurcation 
from wavenumber k = of the linear wave spectrum, shown in Figfl] (b), when there is no 
possible resonance with any linear wave for any value of k. The solitary wave speed v is 
greater than the linear long wave speed, i.e. v > 1, while the speed of a linear wave of any 
wavenumber is smaller, i.e. p < 1. This becomes generally impossible when the symmetry is 
broken. Instead, radiating solitary waves arise for the case when there is a possible resonance 
with the upper branch for some finite non-zero value of k. For example, a possible resonance is 
shown in Figjl] (a) for v = p = 1.3. The solitary wave solutions of Eq. Q, viewed as particular 
solutions of the coupled equations in the symmetric case, constitute a one-parameter family of 
so-called embedded solitary waves (e.g., [67\ I38|). Recently, radiating solitary waves have been 
experimentally observed in two- and three-layered elastic waveguides with soft bonding layers 

From the studies of solitary waves in Refs. [23|. I24|. [25] and our recent studies of the dynamics 
of weakly nonlinear wave packets |40j we know that the dynamics and the asymptotic models 
depend on the relative speeds of the waves. A question arises, to what extent does the difference 
between the characteristic linear speeds of the two wave operators (i.e., c— 1) effect the dynamics 
of the nonlinear waves in Eqs. ([2])? In particular, if we take the initial conditions in the form of 
the solitary waves of the uncoupled Boussinesq equations, will the outcome be different for the 
cases c — 1 = O(e) and c — 1 = 0(1), where e is the natural small parameter of the Boussinesq 
model? The following analysis shows that this difference is crucial. 

3 Weakly nonlinear solution of the initial- value prob- 
lem 

In this section we are concerned with unidirectional waves, which constitute the leading-order 
terms in our asymptotic multiple-scales expansions. We use an averaging procedure which 
allows us to derive asymptotic reductions to simpler asymptotically exact models in the form of 
coupled and uncoupled Ostrovsky equations, and construct a weakly nonlinear solution of the 
initial-value problem in terms of solutions of the derived Ostrovsky equations. 

For these purposes we need to rewrite system (pi) in the original unsealed form (see Ref. [8] ) , 
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substituting 

/ = e/, g = eg, t = yfet, x = 
into the system Q, to obtain (omitting the tildes) 

1 



'€ X, 



e 2 S, 7 



e 2 l 



ftt — fxx — e 



gtt - c g xx 



(f )xx + fttxx - 6(f - g) 



\a{g 2 ) xx + (3g ttxx + 7(7 - g) 



(4) 



As any Boussinesq-type system, system Q appears as an approximation containing 0(1) and 
0(e) terms (see Ref. [5]). 

We consider the Cauchy problem for Eq. Q on the infinite line, imposing the following 
initial conditions: 



f\t=o = F(x), g\ t=0 = G(x), 
ft\t=o = V(x), g t \t=o = W(x). 



(5) 
(6) 



Some local existence results applicable to this problem were recently obtained in Ref. |41j 
(Theorem 2.4 and Remark 2.5, according to Ref. [12]). In this paper, we are concerned with 
the explicit construction of the weakly nonlinear solution of the Cauchy problem in terms of the 
asymptotically exact (KdV-like) models for unidirectional waves. 

We assume that the initial conditions are sufficiently rapidly decaying at both infinities, so 
that to leading order the initial (t = 0(1)) evolution of the Cauchy data is described by the 
classical D'Alembert's solution 



where 



fo{t,x) = f ix-t) + f£(x + t), go(t,x) = g Q ( x - ct) +gv(x + ct), 



g{x ± t) = X - (f{x ±t)± j X±t V{x)dx\ , 

1 / 1 r x±ct \ 

g£ (x ±ct) = - [G(x ±ct)±-J W(x)dx) . 



(7) 
(8) 



In general, / Q and g^ are some step-like functions. In what follows we restrict our considerations 
to the case when these functions are sufficiently rapidly decaying at infinity (i.e. V(x)dx = 
and J™ oo W(x)dx = 0). 

To describe the subsequent (t = O^ 1 )) evolution of the given initial data we introduce the 
slow time T = et and look for the weakly nonlinear solution of the Cauchy problem Q - ^ in 
the form of asymptotic multiple-scales expansions. The form of these expansions depends on 
the difference between the characteristic speeds of the linear wave operators in system Q , and 
below we consider the two main cases, when c — 1 = 0(e) and c — 1 = 0(1). 



3.1 Case I: c — 1 = 0(e) 

In this case, we rewrite system Q as 
ftt fxx — e 



-j{f 2 )xx + fttxx — S(f — g) 



gtt — g% 



1 c 2 - 

2 a (a 2 )xx + figttxx + i(f - g) H - 



(9) 
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where g — L 



~ O(l) since c — 1 = O(e), and look for the solution in the form 

/ = r(£,T) + /+(r ? ,T) + ef l (^r,,T) + 0(e 2 ), 
5 = 5" (£, T) + 5 + (t/, T) + eg\t V, T) + 0(e 2 ). 



(10) 



Here, £ = x — t,rj = x + t,T = et, and we consider each wave in its own reference frame. We 
will view the leading order approximation of the linear solution Q - Q (i.e., for c — 1 = O(e), 
c is replaced with 1) as initial conditions for the functions f~,f + ,g,g + with respect to the 
slow time T, i.e. 

^=0 = ^, 9^=0 = 90 (11) 

(this is later derived at leading order, when we substitute our asymptotic expansions into the 
initial conditions ^ and 

Substituting the expansions (10) into Eqs. ([9]), we find that the equations are satisfied at 
leading order, while at 0(e) we obtain 

- 4 4, 



Mr, 



= (2/f + r/ € - + + (-2# + / + /+ + 

+ 2/-/+ + /+/" + /"/+ - *(/" +f+-g-- g+), 

= (29t + a 9~dJ + + (-29t + a 9 + 9^ + Pa^ m )r, 

+ u(2g-g+ + g + g~^ + g~g} m ) + j(f~ + f + - g~ - g + ) 



(12) 



+ 



c 2 -l 



(% + gi,) 



(13) 



We next average equations (12) and (13) with respect to the fast time t, considering 

i r 

lim - / . . . dt 

r-J-oc r Jq 

at constant £ or 77, i.e. in the reference frame moving with the linear speed of the right- or 
left-propagating waves, respectively. Indeed, requiring that / , g 1 and their derivatives remain 
bounded (which is necessary to avoid the appearance of secular terms in expansions (10)), we 
see that, for example at constant £, 



lim — 

T— >oo T 



fhdt 



lim — 

T— »00 2T 



C+2t 



1 r^n£+ 2r 



lim ^ [/fjf 

t^oo 2r s ? 



0, 



and similarly for g , as well as for f 1 and g 1 at constant 77, showing that the averaging re- 
sults in zeros in the left-hand sides of Eqs. (12) and (13). Similarly, assuming that functions 
/~; / + ) g ^g + an d their derivatives remain bounded and sufficiently rapidly decaying at infinity 
for any fixed T (the assumptions are consistent with the relevant numerical experiments), and 
averaging the entire Eq. (12) with respect to t at constant £, we obtain 





+ lim — 

T^-OO 7~ 



(-2/+ + /+/+ + /+ ), + 2/~/+ + f+fZ + /"/+ - S(f + - g + )]dt 



= (2f T +f-f i+ f^)^-5(f--g-) 

+ lim - J [(-2/+ + /+/+ + + 2/-/+ + + /- /+ - «J(/+ - ff +)]dry 



2/f +r/f -*(/■ 



(14) 



while averaging Eq. ( 13 ) at constant £ results in 



c 2 -l 



0=[2g T +ag g ( +/ 3 % ? J ? +7(/ ~9 H — %• 

Similarly, averaging Eq. ( |12[ ) at constant ?y under the same assumptions yields 

= (-2/+ + /+/+ + f+ m ) v - 5{f + - g+), 



while averaging Eq. ( 13 ) at constant rj results in 



2 1 



= (-20+ + ag + g+ + /?<?+ J + 7 (/+ - 5 + ) + — — g+. 



(15) 



(16) 



(17) 



'v e 

Thus, to leading order we obtain two systems of coupled Ostrovsky equations. Next, substituting 
the Eqs. (14) - (17) back into Eqs. (12) and (13), we obtain equations for the higher-order 
corrections 

fir, = ~\ (2/f # + f + f^ + rfin) , 9\n = -f fasti + <? + % + 9~9^) , 
which imply 

Z 1 = (2/-/ + + /f y /+dr? + /+ y /-c^) + 01 (£, r) + ^(r?, T), 

ff 1 = -| (V<? + + 9{ J 9 + d V + g+ j g'd^j + 2 (£, T) + ^ 2 (r/, T). 

The presence of four arbitrary functions allows us to satisfy not only the equations, but also 
the initial conditions ([5| and ^ up to 0(e 2 ), constructing therefore an accurate asymptotic 
solution of the initial- value problem (within the accuracy of the problem formulation) . 

Indeed, substituting our expansions (10) into the initial conditions §5§ and ^ to leading or- 
der we recover the formulae ( 11 ) for the initial conditions, while at 0(e) we obtain D'Alembert's- 
like formulae for the functions cf)i(£,T) and ipi(rj,T),i = 1, 2: 



MS,?) 



Riiti,T)+ f R i2 (x,T)dx 

J — 00 



Rn(v,T) 



Ri2(x,T)dx 



(18) 



where 



R n {x,T) 
R 12 (x,T) 
R2i{x,T) 

R22(X,T) 



t=0 



+ fx + 7 ( f + fi - rfv + /« / f +dr i - tin / 



t=0 



a 
1 



2g 9 + + % J 9 + dr] + 5+ j g d£ 
9t+9t + J f# + #f ~ 9~9^ + % J 9 + dr] - g+ v J g~d£ 



t=o 



Thus, within the accuracy of the problem formulation (i.e. 0(e 2 )), the dependence of functions 
4> and ip on the characteristic variables is determined, while their dependence on the slow time 
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T is inherited from the dependence of the leading order waves, or it may be neglected, at least 
for sufficiently small values of time. 

The leading order systems of coupled Ostrovsky equations for unidirectional waves can be 
rewritten in a symmetric form in the reference frames moving with the average linear speed 
c = i -e- f° rma lly changing £ and rj in (14) - (17) to £ = £ — AcT and f\ = rj + AcT, 
A 



^2^, which yields 



2(/F -Ac/p + r/f + % 



2(# T + Ac 9? ) + ag g- + fig 



8{f--g-), 
= -7(/~ -5" 



and 



[2(/+ + Ac/±)-/+/+-/+^]_ 
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Systems of coupled KdV equations have appeared in the literature before (see Refs. [251 136] 
and references therein) . To the best of our knowledge this is the first appearance of the coupled 
Ostrovsky equations. 



3.2 Case II: c - 1 = 0(1) 

In this case, we look for the solution in a different form: 



/ = r(£i,T) + / + (??i,T) + ef l (^, m ,T) + 0(e 2 ) 
3 = g-fa,T) + g + (v2,T) + eg 1 ^, t&, T) + 0(e 2 ), 



(19) 



where £1 = x — t,r]i = x + t, and £2 = £ — ct, 772 = £ + ct are the two pairs of characteristic 
variables for the two linear wave operators in Eqs. Q, and again, we consider each wave in its 
own reference frame. 

Substituting our expansions ((19} into Eqs. @ we obtain 



*4 



'/1 



-4c 2 4 



2</2 



= (2/ T + /"./£ + / 6fl€l ) Cl + (-2/+ + /+/+ + /, 

+ 2fgf++f+f^ i+ rf+ iril -6(r + f+-g- 

= {leg? + ag~gg + P^gg^)^ + (-2cgJ + a5 + s>+ + /3c 2 #. 

+ a(2%9^ 2 + 9 + % 6 + 9~9^ 2m ) + 7(/~ + f + ~ 9~ ~ 9 + )- 



(20) 



(21) 



Under the same assumptions as before, we can average Eqs. (20) and (21) with respect to 
t, considering 

1 r 

lim — / . . . dt 

T^OO T Jq 

at constant £1 or 7/1, and £2 or 772, respectively, obtaining in this case four uncoupled Ostrovsky 
equations: 



(2/ T — f + f^ — frjirjirji)m = 



(2^ - ac/ + 94 - Pc 2 g. 



V2V2V2' r l2 



-19 



(22) 
(23) 
(24) 
(25) 



S 



and equations for the higher-order corrections 

4m = ~\ ( 2 €/m + f + f[^ + f'fmm) ~ |GT + 9 + ), 
dim = "^2 ( 2 %4 +9 + 9^ 2 +9'9 mm ) ~ + / + ), 

where in the right-hand sides we have solutions of the leading order Ostrovsky equations 
(1 + c)Ci + (1 - c) m 



(26) 
(27) 



g-fa,T)=g 
and 

r(£i,r) = / 



+ , r . I , (1 - C)gl + (1 + C) m ,„ 

9 {V2,T) = I ,T 



(c + l)6 + (c- l)^ 
2c 



Remarkably, the particular solutions of Eqs. (26) and (27) are bounded functions, because 



of the special property of the smooth solutions of the Ostrovsky equation, namely 



f-dti = 0, 



f+d Vl = 0, 



g d^2 = o, 



g + dt]2 = o. 



(28) 



Indeed, the solution of Eqs. (26) and (27) can be found in the form 



/ 1 = -|(2r/+ + / 6 J /+<% + /+ f r^i) +f P i(&,T) + f p2 (r ]2 ,T) + MZi,T) + MvuT), 

r 9 + + g^J g + dm + g^2 J 9~d^ +g P i(£i,T) + g P 2(vi,T) + ^2(6,^) + ifj 2 (v2,T), 



9 l = ~ 
where 



4? 



/i 



pi 



9 P i 



6 r v 



g (u,T)dudv, / 3 



p2 



—00 j —00 



1 



/ (u,T)dudv, g p2 



—00 t/ —00 
m rv 



00 J —00 



1 



g + (u, T)dudv, 
f + (u, T)dudv. 



00 *j —00 



Let us consider f p \, for example. Here, g~{u, T)du = 0, because of the mentioned property 
of the Ostrovsky equation, immediately obtained by integrating Eq. (24). Moreover, using Eq. 
(24) and recalling that solutions are decaying at infinity, we obtain 

1 f°° 

- / (2cc/t + ag-g- + (5c 2 g~ vv )dv 

7 J — CO 



oo rv 



g (u,T)dudv 



OO J — oo 



7 

2c d 
"7 cff 



<? (w, T)dv = 0. 



Therefore, the particular solution f p \ is a bounded function, and limg 2 _ > -|- 00 f p \ = 0. Similarly, 
other particular solutions are also bounded functions and there are no secular terms. 

The presence of four arbitrary functions allows us to satisfy the initial conditions with the 
desired accuracy, as in the previous case. Substituting our expansions (19) into the initial 
conditions ^ and ^ to leading order we again recover formulae (11), while at O(e) we obtain 
D'Alembert's-like formulae for c/>j(£j,T) and ^i(r]i,T),i = 1,2: 



Rii(£i,T) + 



Rn(r]i,T) 



1 /* 



i r 



R i2 (x,T)dx 



Ri 2 (x,T)dx 



(29) 



9 



where c% = 1, c% = c and 



Ru(x,T) 



Ru(x,T) 



R22(X,T) 



t=0 



„2 _ 1 



x fV 



[g (it, T) + g + (u, T)] dudv, 



oo J — oo 



y [5 K T ) ~9 + {u,T)] du, 



t=o 



a 

4<? 



+ 



7 



25 5 + % y 5 ^2 + g m J g d& 

t r r [r(«,T)+/ + (u,r)]d«di;, 

J —00 J — 00 

- T [/-(«, r) -/+(«, 

J —oo 



t=0 
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Thus, in both cases c— 1 = 0(e) and c — 1 = 0(1), the asymptotic multiple-scales expansions 
and the averaging procedure described above have allowed us to construct nonsecular weakly 
nonlinear solutions of the given initial-value problem for the values of time up to 0(e _1 ), within 
the accuracy of the problem formulation. To construct a more accurate solution, and for greater 
values of time, one would need to know higher-order terms in the orig inal cRB Eqs. Q. 
However, the derived hierarchy of Ostrovsky equations will still describe leading order terms in 
these expansions, making the study of the long-time evolution of its solutions interesting. 

To finish this section, we would like to make an important comment that although smooth 
solutions of the Ostrovsky equation must satisfy the zero mass constraints ( |28| ) (and similar 
conditions /^(Z - — g~) d£ = 0, etc. in the case of coupled Ostrovsky equations of section 
3.1), this does not impose any forbidding restrictions on the choice of the initial conditions §5§ 
and ^ for the cRB system Q. Indeed, initial conditions for the Ostrovsky equation can always 
be modified by adding a long but very small amplitude (i.e. 0(e 2 ) or smaller) 'pedestal' (e.g., 
0(e n ) constant over the finite 0{e~ n ) interval, n > 2), so that the composite initial condition 
has zero mass, but this does not lower the accuracy of the asymptotic solution; and a transition 
to the zero mass solution is very fast (see the discussion in Ref. |43j ) . Numerical simulations for 
the Ostrovsky equation show that if a smooth initial condition has a nonzero mass, numerically 
the solution adjusts immediately, since this initial condition can be viewed as an approximation 
to the composite zero mass solution (with any given accuracy). Note that similar issues appear 
in connection with several other equations, for example, the Kadomtsev-Petviashvili equation 
(see Ref. [S]) and, more recently, the short-pulse equation (see Ref. [45]), where the notion of 
the initial time layer has been introduced to describe such transitions [H] (see also Ref. [13] 
and references therein). 



4 Numerical simulations 

In this section we discuss numerical simulations of solutions of system Q. We implement a 
finite difference scheme which is an extension of the scheme developed in Ref. [36] for a single 
regularised Boussinesq equation. Our emphasis in this section is to compare the numerical 
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solutions for the two cases discussed in the previous section, i.e. when the difference in the 
characteristic speeds of the system is of O(l) or 0{e). 

We let x G [—L,L], for finite L, and discretise the (x,t) domain into a grid with spacings 
Ax = h and At = k. The solutions f(x,t) and g(x,t) of Eqs. ([2J are approximated by the 
solution f(ih,jk) and g(ih,jk) (for i = 0,1,..., iV and j = 0,1,...) of the difference scheme, 
denoted fij and g^j. 

Substituting central difference approximations into system ^ we derive the following dif- 
ference schemes for the system ([2]): 

- /i-ij+l + (2 + h 2 )fij + i - /i+ij+i = (k 2 - 2)[/i_ij - 2 fij + fi+ij] + 2h 2 fij 

+y[(/«-y) a -2(/y) 2 + (W 2 ] 

— (2 + h 2 )fi t j-i + /j+ij-i 
-h 2 k 2 S(f id - g itj ), (30) 



+ (2/3 + h 2 )g iJ+1 - (3g i+ i d+1 



(k 2 c 2 
ak 2 



W)[9i-l,3 ~ 2 9i,j + ffi+ij] + 2/i 2 5i,i 



+^fl , i-lj-i - 



(2/3 + /i 2 )5i,i-i + Pgi+ij-i 

-9i,i)- (31) 
The boundary conditions are imposed far enough from the propagating waves, thus we set 

fo,j = Inj = goj = gNj = 0, Vj. (32) 

The initial conditions are chosen in the form of the co-propagating pure solitary wave solutions 
of the uncoupled equations (5 = 7 = 0): 



fi,0 



9i,o 



Ai sech 



Ai sech 



x 

X 

A2 



fi a = Ai sech 



gi,l = A 2 sech 



x — v\k 

x — v 2 k 
A9. 



Vi, 



(33) 



where A\ = 3(vf - 1), A 2 



*(v 2 



Ai 



l)-3 and A 2 = 2 V2 V?(f 2 2 - c 2 )" 



The nine point implicit difference schemes (30) and (31), with tri-diagonal matrices of constant 
coefficients, are solved simultaneously using a Thomas Algorithm (e.g., |47j ) . 

A single Boussinesq equation with arbitrary coefficients (i.e. system ([2]) for 5 = 7 = and 
/ = 0) is used to examine the scheme's stability. This solution can be approximated by the 
solution of the difference scheme (31) with 7 = 0. We first linearise this scheme by setting 
9i,j =5o + 9i,j where go is a constant such that go > gij Vi, j. Using a Von-Neumann stability 
analysis we substitute gij = G^e l9th (where i 2 = —1) into the linearised version of (31) with 
7 = and derive 



G 2 - 2fj,G + 1 = where p = 1 



2k 2 (c 2 + a#o)sin 2 
/i 2 + 4/3sin 2 f 



(34) 



For stability we require \G\ < 1V9 and arbitrary k, h, which is true provided |/x| < 1 and thus 
implies k < k c = \J ^2^g - Hence the roots of the quadratic in G have modulus one and the 
linearised form of the difference sheme (31) with 7 = is stable provided k < k c . (In practice, 
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we used a stricter condition k < ^k c , to accommodate for the effects of nonlinearity) . It can be 
shown that the principal truncation error of this scheme is 0(h 2 k 2 ). 

Numerical simulations for the symmetric case (scheme (30) with 5 = 0) compared with the 
known analytical solution Q reveal that the choice of discretisation can reduce the maximum 
error, across x for a given time, to as low as O(10 -5 ). This accuracy is within the range 
which is deemed suitable from previous work on Boussinesq-type equations (see Refs. |48j - 
|54j ) . The step size h = k = 0.01 results in errors of this order and is thus chosen for our 
simulations. We also numerically approximate u(x,t) and w(x,t) via Simpson's rule using the 
relationship J_ L f,gdx = u, w, due to the boundary conditions (32 ), and hence utilise the energy 
conservation law given in Section 1, although integration introduces additional errors. The 
conserved quantity A2dx was monitored and for simulations with h = k = 0.1 the energy 
was conserved, within the chosen time interval, up to 0.021% and 0.006% for the results shown 
in Fig. 2 and Fig. 3, respectively. For smaller step sizes these computations become very time- 
consuming, but there are no noticeable differences in the plots of the solutions for h = k = 0.1 
and h = k = 0.01. (A useful discussion of the difficulties associated with the accuracy of 
conservation laws in finite-difference schemes can be found in Ref. |19|.) 
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ure 2: Generation of a radiating solitary wave for c = 1.05, a = /3 = l,7 = <5 = 0.01; v± = i>2 = 1.3 fr 
solitary waves of the uncoupled equations. 



The emergence of radiating solitary waves replacing the initial pure solitary waves in both 
components of / and g is shown in Fig(2] for the case c = 1.05, f3 = a = 1,5 = 7 = 0.01; v i = 
i>2 = 1.3. This result agrees with the discussion in Section 2 and numerical studies in Ref. [8]. 
The numerically determined wavelength of the oscillatory tail at t = 300 is (36.5 ± 0.1) for / 
and (37.0 ± 0.2) for g, which is close to the theoretical prediction (ss 36.7) for p = 1.3 using the 
dispersion relation (see Fig.l). 

For the case c = 2, /3 = a = 1,5 = 7 = 0.01; v\ = 1.3, = 2.3, the initial solitary waves 
are replaced by dominant wave packets in both components, shown in Fig(3} The emergence 
of a wave packet in both / and g can be observed almost instantaneously and for long time 
appears stable. From our asymptotic analysis for c — 1 = 0(1) in Section 3 we found that to 
leading order the solution for both / and g, for right propagating waves, are the solution to 
the two Ostrovsky equations (22) and (24). At t = 600 in Fig J3] the leading wave packet in g 
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Figure 3: Generation of wave packets for c = 2, a = /3 = 1, 7 = S = 0.01; v\ = 1.3, ^2 = 2.3 from pure 
solitary waves of the uncoupled equations. 



is qualitatively very similar to the numerical solution of the Ostrovsky equation studied in Ref. 
[5]. Similarly, at t = 600 for / a similar but smaller and slower moving wave packet is present 
around x = 600. In Ref. [5] it is shown that a parameter 00, equivalent to ajo = 6 ^/^ an d 

a g o = for / and g respectively in our system, determines the emergence of this distinct 

wave packet. For our simulations in Figj3j this parameter lies within the range for which the 
distinctive wave packet will emerge and hence it can be seen in both / and g. As afo,a g o —> 
there exists a range for which there is no emergence of a distinctive wave packet. Alternatively 
as the parameters are increased, the faster the wave packet emerges and in the case of g, the 
faster it will move away from the rest of the solution. The range for our results for which this 
transition occurs is also in very good agreement with the results for the Ostrovsky equation in 
Ref. 0. 

The simulations in this section confirm that there is a distinctive difference in the qualitative 
behaviour of the solution upon varying the difference in the characteristic speeds, and hence 
support our asymptotic analysis. From our numerics we can conclude that for pure solitary 
wave initial data in the parameter range c — 1 = O(e), stable radiating solitary waves emerge in 
both components / and g. However as we increase the difference in the characteristic speeds, 
i.e. let c— 1 = O(l), we see the emergence of wave packets, agreeing with the numerical solution 
of the Ostrovsky equation [5J, in both components / and g. 

5 Conclusions 

In this paper we addressed the question of constructing a weakly nonlinear solution of the initial- 
value problem for coupled Boussinesq-type equations for localised or sufficiently rapidly decaying 
initial data, generating sufficiently rapidly decaying right- and left-propagating waves. Crucially, 
we considered the general case, when the two linear wave operators have different characteristic 
speeds, which complicates the analysis since in this case the number of characteristic variables 
(four) is greater than the number of independent variables (two). Further generalisations to the 
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case of more than two equations (and characteristic speeds) are straightforward. 

We introduced two different types of asymptotic multiple-scales expansions for the cases 
c — 1 = 0(e) and c — 1 = 0(1) and averaged with respect to the fast time, which allowed us to 
derive to leading order a hierarchy of asymptotically exact coupled and uncoupled Ostrovsky 
equations for unidirectional waves. We then constructed the nonsecular solution of the initial- 
value problem in terms of solutions of the derived leading order models for the values of time 
up to 0(e _1 ), within the accuracy of the problem formulation. To construct a more accurate 
solution, and for greater values of time, one needs to know higher-order terms in the original 
cRB Eqs. Q. 

We performed numerical simulations of the original unapproximated coupled regularised 
Boussinesq Eqs. Q for the initial conditions in the form of co-propagating pure solitary waves of 
the uncoupled equations, and compared our numerical results with the known numerical results 
for the Ostrovsky equation [5], which confirmed predictions of our leading order asymptotic 
theory. Expanded numerical studies for this and other types of initial conditions are currently 
underway. 

The approach developed in this paper is generic and can be used to construct weakly nonlin- 
ear solutions of some other initial-value problems, and in other physical contexts. In particular, 
it is interesting to derive a hierarchy of Ostrovsky equations and consider the initial- value prob- 
lem for the original area of appearance of this equation (rotating ocean) , which will be discussed 
somewhere else. 

When (5 = 7 = 0, our solution yields an explicit weakly nonlinear solution of the initial-value 
problem for the single Boussinesq equation 



fit fx 



(f 2 U + f 



ttxx 



f\ 



t=0 



F(x) 



ft 



t\t=o 



V(x), 



(or any other asymptotically equivalent form of this equation) for the case when the ini- 
tial conditions generate sufficiently rapidly decaying right- and left-propagating waves (i.e. 
J^oo V(x)dx = 0). The solution has the following form: 



f = f-(ti,T)+f+( V ,T)+e 



\(^f~f + + f[ j f + d7] + f+ j rdd)+^T)+^T) 



where £ = x — t,rj = x + t,T = et. The functions / and / + are solutions of the initial-value 
problems for two Korteweg-de Vries equations 



1 



1 



1 



1 



0, 



/ |t=o 



/x±t 
V(x)d 
-oo 



integrable by the Inverse Scattering Transform [55] (see also Refs. [12\ [56] and the relevant 
discussion of matching of the near-field and far-field solutions for unidirectional waves in Ref. 
|57] ) . The functions </> and tp are given by the formulae 



</>(£, T) 
where 



P(£,T)+ f Q(x,T)dx 

J — oo 



P(rj,T) 



Q(x, T)dx 



P(x,T) 



Q(x,T) 



2f-f + + f7 / f + dr, + f+ / f-di 



t=o 

Jr + It + ~ [ r.i'7 r/,7 + / f + dv ~ / J~< 
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(within the accuracy 0(e 2 ) of the problem formulation, the dependence of the functions 4> and 
ip on T is inherited from the dependence of the leading order functions /~ and / + , or it may 
be neglected, at least for sufficiently small values of time). 

For the practical applications of the constructed solution it is useful to remember that within 
the accuracy of the problem formulation in Q (i.e. 0(e 2 )), the initial conditions ^ and ^ 
can be represented in the form 

f\t=o = F°(x) + eF\x) + 0(e 2 ), g\ t=0 = G°(x) + eG\x) + 0(e 2 ), 
ft\t=o = V°(x) + eV\x) + 0(e 2 ), g t \t=o = W°(x) + eW\x) + 0(e 2 ), 

which not only allows one to formally satisfy the zero mass constraints for , by adding 
appropriate 0(e 2 ) 'pedestal' terms, as explained in section 3, but also gives us some flexibility 
with the choice of initial conditions for the auxiliary IVP problems for unidirectional waves by 
splitting the functions F(x),G(x),V(x),W(x) into a 'nice' O(l) part (i.e. such that the IVP 
problems have some favourable analytical properties, e.g. from the viewpoint of the Inverse 
Scattering Transform when these are for the KdV equations) and a small O(e) remainder, which 



can be readily accounted for in D'Alembert's-like formulae (18) and ([29j) for the functions 4>i 
and tpi, i = 1,2. 

Finally, we would like to emphasise the importance of the Ostrovsky equation as a canonical 
asymptotically exact model, similar to the Korteweg-de Vries model. The reduced form of the 
Ostrovsky equation 

(r]t + vrjVx)x = At? 

was recently shown to be an integrable equation [58, 59J, reducible to the Tzitzeica equation [60] , 
We believe that the full Ostrovsky equation also might have some 'nice' analytical properties 
(although it is not necessarily integrable in the conventional sense). 
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